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ABSTRACT 

In order to explain the main characteristics of the observed population of ex- 
trasolar planets and the giant planets in the Solar System, we need to get a clear 
understanding of which are the initial conditions that allowed their formation. To this 
end we develop a semi-analytical model for computing planetary systems formation 
based on the core instability model for the gas accretion of the embryos and the oli- 
garchic growth regime for the accretion of the solid cores. With this model we explore 
not only different initial discs profiles motivated by similarity solutions for viscous 
accretion discs, but we also consider different initial conditions to generate a variety 
of planetary systems assuming a large range of discs masses and sizes according to 
the last results in protoplanetary discs observations. We form a large population of 
planetary systems in order to explore the effects in the formation of assuming different 
discs and also the effects of type I and II regimes of planetary migration, which were 
found to play fundamental role in reproducing the distribution of observed exoplanets. 
Our results show that the observed population of exoplanets and the giant planets in 
the Solar System are well represented when considering a surface density profile with 
a power law in the inner part characterised by an exponent of —1, which represents a 
softer profile when compared with the case most similar to the MMSN model case. 
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1 INTRODUCTION 

So far 464 planets have been found orbiting stars in the solar 
neighbourhood. This population of planets, though strongly 
biased to those planets that are easier to be detected with ra- 
dial velocities, present tendencies that the theoretical mod- 
els of planetary formation must reproduce. With the end of 
explaining the main characteristics of the observed popula- 
tion of extrasolar planets, including the planets in the Solar 
System, we need to get a clear understanding of which are 
the initial conditions that allow their formation. 

The fact that the planetary systems formed from a disc- 
like nebula was recognized centuries ago, but even today our 
understanding of this processes and initial conditions remain 
uncertain. With this question in mind, the main objective 
of this work is to determine which is the nebula model that 
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allows us to form planetary systems whose planets reproduce 
the observational sample of extrasolar planets. 

In recent years there have been s everal sem i ana- 
lytical models of planetary formation dlda fc Linl 12004 
iMiguel fc Bruninil 120081 ; iMordasini et al.1 l2009h and plane- 



tary systems formation (|Miguel fc Bruninil l2009'). which in- 
tend to reproduce the observational sample of exoplanets 
and explain with a simple model the main characteristics of 
this distribution in order to get a better understanding of 
the process of planetary formation. For simplicity the sur- 
face density model assumed in these models as basic build- 
ing blocks for the planetary formation process are based on 
conventional power-law models, as t he standard m inimum 
mass solar nebula (MMSN) model of lHavashil i|l98ll ). which 
is an approximation to the gas-dominated nebula, taking 
into account the primitive composition of the solar disc. This 
model, although efficient on its simplicity, suffers from mul- 
tiple disadvantages. 

The minimum nebula fedgeworthlll949l; lKuip"erlll956l ; 
ISafronovl [l969l : IWeidenschillind Il977l : iHavashil 119811 ), was 

constructed by assuming that there should be added to the 
masses of the present planetary bodies enough icy materi- 
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als and hydrogen and helium to achieve solar composition 
and then this mass should be smeared over the surround- 
ing place. So the model was based on the strong implicit 
assumption that planets were formed at their present loca- 
tions, they accreted all of the solids in their vicinity and the 
accretion was perfect, which means that the formation of one 
giant planet did not have consequences on the formation of 
the others. Despite this model is a reasonable good approxi- 
mation in the intermediate region of the solar nebula, posses 
too much surface density in the inner region, under predicts 
it in the outer nebula, and the total mass of the disc is infi- 
nite in this formulation, so inner and outer boundaries must 
be specified. For all these reasons and despite its widespread 
use, it is clear that this model must be updated or changed. 

With the spirit of retaining a simple model for the sur- 
face density in the nebula and get a more realistic model 
consistent with protoplanetary discs observations, we con- 
sider a nebula with a surface density profile motivated by 
si milarity solutions for v is cous accretion discs as shown 



Lynden-Bell fc Pringld i|l974 ) and in lHartmann et al.l 



(1998). This simple formula, which is characterised by a 
power-law in the inner part of the disc and an exponen- 
tial decay in the outer part, is a reasonable alternative for 
use in protoplanetary discs because it is simple, it matches 
the protoplanetar y discs observations l|Andrews et al.|[2009l ; 
llsella et al.l 120091 ) and unlike power-law approximations, 
these surface density model predicts a sharp outer edge to 
the nebula without introducing an arbitrarily cut to the disc. 
There are several free parameters in the nebula model that 
we have to adjust with the aim of explaining the distribu- 
tion of observed exoplanets and the giant planets in our 
Solar System, such as the exponent of the power-law in the 
inner part of the disc, whose different values will bring a 
significant impact on the formation of giant planets as we 
discuss in this work. 

So with the aim of finding the planetary system for- 
mation model that best fits the observations (including 
our Solar System), we have developed a semi analytical 
model where the main d ifference with our previous model 
jMiguel fc Bruninil 120091 ) is the initial nebula profile. The 
model, which is shown in section [21 is based on the nucle- 
ated instability model, where the solid cores grow in the oli- 
garchic growth regime. We also include the effects of having 
embryos in a gaseous disc, considering type 1 and II regimes 
of planetary migration, which were found to be fundamen- 
tal in reproducing the distribution of observed exoplanets. 
The results are shown in section [3] where we also show the 
comparison with the observations and the summary and con- 
clusions are shown in section [4] 



2 MODEL AND BASIC EQUATIONS 

In this section we show the model co nsidered, which is 
based on a model previously developed bv lMiguel fc Bruninil 
l|2008l . l2009l . l2010l i. As the main goal of this work is to show 
which are the consequences of considering different density 
profiles and which one suits better the observations and 
the Solar System, we will explain the protoplanetary nebula 
model in detail. 



2.1 Protoplanetary Nebula Structure 

The structure of the protoplanetary nebula is usually as- 
sumed to follow a power-law distribution in the form, 
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where Eo is the surface density at the arbitrary radius do - In 
this formulation, inner and outer boundaries of the disc must 
be arbitrarily specified otherwise the total mass of the disc 
would be infinite. This is one of the limitations suffered by 
this formulation, that was initially motivated by models as 
th e minimum ma ss solar nebula model (MMSN) developed 
bv lHavashil (|l98lh . 

In order to consider a more app ropriate disc structure 
mo del, we are based o n the works of lAndrews et all (|2009h 
and llsella et all |2009l ). who adopted a density profile char- 
acterised by a power-law in the inner part of the disc and an 
exponential decay in the outer parts of the disc and is based 
on the similarity solutions of the surface density of a thin 
Keplerian di sc subject to the gravity of a. point mass (M») 
centr al star (|Lvnden-Bell fc Pringlel [l974l ; lHartmann et all 
1998). In this formulation the gas surface density is given 
by, 



E a (a) = E° ± e 
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where a c is a characteristic radius beyond which the density 
drops sharply, 7 is the exponent that defines the density 
profile and E° is a parameter that is determined from the 
total mass of the disc which is given by, 
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solving the integral we get the expression for Eg, 



(4) 



We adopted a temperature profile for a disc optically 
thin, given by 



(5) 



with Lq and the Sun and Stellar luminosity. The volu- 
metric density of gas is 

_(z) 

p = p e w (6) 

where h is the disc scale height (h ~ 0.05(a/lait) 1//4 a) and 
po is the volumetric density of gas in the midplane of the 
disc given by 



po = 8.33 x 10" 



, , \ !/2 / \ -7-5/4 



/ a c 



.(^) 2 -Vl \ 1/8 _ 3 

LI gCm 



(7) 



We will assume that the solids surface density has an 
expression similar to equation [2] 
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Figure 1. Solid surface density for three discs with the same total 
mass, characteristic radius and metellicity but different values for 
7. The dotted line represents 7 = 1.5, the dashed line is for 7 = 1 
and the solid line represents the density for 7 = 0.5. 



where r/i ce is a function introduced in order to represent 
the change in the solids beyond the radius where the water 
condenses, 
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where the snow line is located at a,i ce — 2.7(-^^) au from a 
central star of mass M* and if zo is the primordial abundance 
of heavy elements in the Sun and we assume the same for 
the disc, then 



= z 



(10) 



following the result o f lLoddersI (|2003l ) we assumed that zq 
0.0149 and for a star of metallicity [Fe/H] then 
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The metallicities are taken random considering that they 
follow a log-normal distribution fit ted from the resu l ts of t he 
CORALIE sample, as is shown in iMordasini etafl i|2009l ). 

In order to illustrate we show in Figure [T]the solids sur- 
face density as a function of semi major axis for three discs 
with a c = 50 au, total mass of 0.05 Mq, solar metallicity 
and different values for 7. The dotted line represents the 
solid surface density for 7 = 1.5, the dashed line shows the 
curve for the case when 7=1 and finally the solid line is 
the solid surface density for 7 = 0.5. We note that as the 
value of the exponent 7 decreases, the amount of solids in 
the inner part of the disc is reduced, while it increases in the 
outer part, fact that brings deep consequences on the forma- 
tion of giant planets, as we will see in the results. The case 
most similar to the MMSN model is when 7 = 1.5 where we 
note an excess of solids in the inner part of the disc, which 
drops sharply beyond the snow line. This case is less favor- 
able for the formation of several giant planets in the same 
planetary system, since the mass is concentrated and as a 
consequence the formation of a single giant planet in the 
snow line is allowed leaving no residual material for others 



after i ts formation. Recently iGuilera. Brunini fc Benvenutol 
have shown that when employing power-law discs 
E oc a~ p , the value p = 1.5 leads to a quickly formation of 
Jupiter that could inhibit the formation of Saturn when the 
formation of both is considered simultaneously. They also 
showed that smoother surface density profiles (p=l, 0.5) 
favoured the simultaneous formation of Jupiter and Saturn. 

The disc is extended between a;„t and a ex t, where the 
inn er radius i s calc ulated accorded to an expression given 
by IVinkovicl (|2006t ) who found it through observations in 
protoplanetary discs, 
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with T 3U b the dust sublimation temperature taken as 
1500°^ and L* and Lq are the stellar and Sun luminosity 
respectively. The inner radius takes values around 0.01 au. 
The outer radius is the radius that contains 95% of the disc 
mass (the total disc mass is aproximately the gaseous mass 
of the disc), 
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where we assume that the total disc mass is aproximately 
the gaseous mass of the disc. Then, 
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which for a disc characterised by 7 = 1 ranges between ~ 90 
and ~ 300 au. 

The total mass of the discs and their characteristic ra- 
dius are taken random considering that they also follow a 
log- Gau ssian distribution tha t we fitted follow i ng th e re- 
sults of lAndrews et"aH l|2009h and llsella et all (|2009T ). As 
the masses of the generated discs could be large enough to 
undergo in a gravitational instability we check that our discs 
are stable. The gravitational stability of a Keplerian accre- 
tion disc wit h sound speed c s is measured by the Toomre 
Q-parameter l|Toomrelll964| ). which is defined by, 

n C a Qk 
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In the case of our disc model, 
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where a value of Q ^ 1 represents an unstable disc. 

As we note in the previous equation, the parameter 
of instability depends on the semimajor axis, therefore it 

changes through the disc. The minimun of equation [16] is 

,7 . 1 

found when a — ( | - j 2-7 a c , so we chose Qa min (the value 

of Q when representative value for the disc. 

Those discs with Qa min > 1 will be stable all over the disc. 
We also notice that the gravitational stability parameter de- 
pends on 7, hence different values of 7 lead to different values 
of the Toomre parameter and therefore two discs with equal 
mass and characteristic radius could be stable or unstable, 
depending on the density profile. 

This is shown in figures [2] where we show the mass vs. 
the a c of all the stables discs generated, and the different 
figures show the results for three different values of the pa- 



rameter 7. Figure 2(a) show the stable discs generated with 
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Figure 2. Mass and characteristic radius of all the discs gen- 
erated with Qamin > 1- the first figure shows the results when 
we consider 7 = 1.5, figure [2(b) | shows the stable discs generated 
when 7 = 1 and finally stable discs with 7 = 0.5 are shown in 
figure Egg] 



a density profile corresponding to 7 = 1.5, 2(b) are the discs 
with 7 = 1 and 2(c) when 7 = 0.5. We note that for larger 



values of 7 we need more mass to obtain a gravitationally 
unstable disc. This suggest that, as is seen in equation 1151 
the global nature of the instability depends both on the mass 
of the disc and on how this mass is distributed. 



Finally, in figures 2(a) and 2(b) we note the presence of 
extremely massive discs. These discs should not be consid- 
ered as keplerian and equation 1151 does not apply to them. 
In order to avoid these very massive discs, we check that the 



disc mass is less than 20% of the mass of the central star 
ijHartmann et al]|l998l ). Thus we consider discs with masses 
up to 0.28 M . 

In order to represent in a simplified way the depletion 
of the gaseous disc, we assume that this mechanism can be 
modeled considering an exponential decay for the mass of 
the gaseous disc, which occurs on time scales between 10 6 
and 1 7 years in agreement to observation of circumstellar 
discs (|Haisch et al.ll200ll ; lHm"enbrandll2005r ). 



2.2 Growth of the Protoplanetary Embryos 

Once we have the protoplanetary nebula defined, we locate 
the initial cores. The first one is located at a,i nt and the 
others are located at 10 Hill radii away from each other until 
it reaches the disc outer radius a ext . Each initial embryo 
has the minimu m mass necessary for starting the oligarchi c 
growth regime ijKokubo fc "idal Il998l : llda fc Makinol ll993T l , 
which is the growth stage where our embryos grow. 
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with m the effective planetesimal mass. 



2.2.1 Solid Accretion into the cores 

The cores grow due to the accretion of solids, gas and also 
due to the collision with other embryos. The rate at which a 
core accu mulates solids in the oligarchic growth regime was 
found bv lSafronovl (1969) and has the form, 
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where f2 is the Kepler frequency, R p and M t are the planet's 
radius and total mass (solid + gas) and a is the velocity dis- 
persion which depends on the eccentricity of the planetesi- 
mals in the disc. We assume that the rms eccentricity of the 
planetesimals in the disc are damped and have reached an 
equilibrium value which is, 

m 101/5^5^1/3^/5 W 

as is shown in lThommes et all (|2003l ) . where p m is the plan- 
etesimal bulk density. With this expression the solid accre- 
tion rate including the evolution of the planetesimal rms e 
and i has the form, 
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where pu is the embryo bulk density, which is equal to the 
planetesimals density, pu = pm = p. 

The solids accretion ends when E s is zero in their feed- 
ing zones. The solids surface density diminishes due to a 
combination of two factors: 

• the cores ate the solids in their feeding zones, 

• the ejection of planet esimals diminished the solids sur- 
face density in the region (jThommes et al 1 l2003l ; llda fc Linl 
12004 ). 
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We also consider that when a core is able to retain gas, 
the drag effect caused by the gaseous envelope on the plan- 
etesimals increases the colli sion cross section of the embryo. 
Following [Chambers! < |2006h . we assume that the enhanced 
collision radius (RcoMsion) of the embryo is given by, 
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where c is the velocity of light, P is the orbital period of 
the protoplanet, n is the opacity of the atmosphere which 
is considered as ~ 4cm 2 <? _1 , fJ. — 2.8 is the mean molecular 
weight of the atmosphere (assumed to have a solar composi- 
tion) and the embryos equilibrium eccentricity is considered 
as ~ 2 in this expression. 

As was said in the beginning of this section, colli- 
sions represent an important evolutionary process which 
plays a significant role in determining the final mass and 
spin state of the planets. The model and consequences of 
embryo-embryo collision s was explained in a previous work 
([Miguel fc Brunimll201(t l. here we say that the result of the 
collisions are the perfect accretion of the embryos involved. 



2.2.2 Gas accretion onto the Cores 

The cores have an associate envelope if the molecular 
velocity is smaller than the escape one. According to 
Mizuno et all (|l978h . once a protoplanet becomes greater 
than ~ the Moon's mass, the core attracts the neighbor- 
ing gas and an envelope forms surrounding it. In the early 
stages of giant planets formation, the gravity is balanced by 
the pressure gradient which is maintained b y the potentia l 
energy released by incoming planetesimals (Mizuno 1980), 
so the stability of this envelope depends on the mass of the 
protoplanet and on the solids accretion rate. When the mass 
of the embryo becomes greater than a critical value, or there 
are no more solids available on its feeding zone, the envelope 
can no l onger be in hydrostatic equilibrium and beg ins to 
collapse (|Mizundll980l ; lBodenheimer fc Pollacklll986l i as a 
result, the gas accretion process began. T he core' s critical 
mass (|Stevensonlll982l ; llkoma et al.ll2000h is given by, 



M cri t ~ 10 
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(22) 



IQ-^M^yr- 1 

iFortier et all (|2007f ) and iFortier et all (120091 ') studied 
the formation and evolution of a protoplanet in situ, in 
the frame of the nucleated instability model. They im- 
proved a numerical code b ased on a Henyey technique 
([Benvenuto fc Bruninil [2005;), and consider that the oli- 
garchic growth regime is the stage for the accretion of the 
solid cores and considered the gas drag effect acting on 
planetesi mals inside th e plan et atmosphere. We used the 
results of IFortier et al.l |2007) introducing an analytic ap- 
proximation to their nu merical results as is explained in 
iMiguel fc Bruninl ([2008), and update here the gas accre - 
tion rate according to their last results (|Fortier et al.l l2009). 
Then the gas accretion process occurs on a rate, 
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where M g is the mass of the surrounding envelope and r g is 
its characteristic growth time, 



r. = 8.36x10* 
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We assume the next limits for ending the gas accretion 
process: 

• the growth of the envelope ends when the planet con- 
sumes all the gas available on its feeding zone, 

• or it opens up a gap on its orbit, 

• and finally the gas accretion ends when -4^ > ' 



lQOyears 



2.3 Planetary Migration 

The gravitational interaction between the embryos and the 
gaseous protoplanetary disc leads to an angular momen- 
tum exchange between them and as a consequence an or- 
bital motion or migration of the embedded embryos occur 
(|Goldreich fc Tremainelll98Cl ). For low mass planets the an- 
gular momentum flux injected into the disc is negligible 
when compared to the viscous transport of angular momen- 
tum, i n this case we have the type I migration regime, fol- 
lowing ^na^eTaD (|2002h the migration rate is, 



migl 



c migl [2.7+1.10 ( jf- 
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as the time scale for type I migration is inversely propor- 
tional to the mass of the disc and the planet, it can be much 
shorter than the disc lifetime, so the factor c mig i is intro- 
duced for considering effects that might slow down or even 
stop migration without introducing a mayor degree of com- 
plexity to the model and /3 has the form, 



dlog(T, g ) 
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When the planet reaches the mass necessary to open 
up a gap on its orbit, the angular momentum flux from the 
planet locally dominates the viscous flux and a new regime 
of planetary migration begins. This is the type II regime , 
characterised by the n ext migration rate l[Lin et all [l99rj ; 
iLin fc Papaloizoulll985l ). 
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with a = 10 3 a dimensionless parameter which charac- 
terises the viscosity, R m = 10e T <«sc au and Tdi SC the disc 
depletion time-scale. 

We assume that both migration mechanisms stop when 
the core reaches the inner edge of the disc. 



3 RESULTS 

We perform a series of simulations varying the density profile 
and the type I migration rate and analyze which are the 
consequences on the planetary formation of considering a 
different prescription for the protoplanetary nebula profile 
and migration rates. In each simulation we generate 1000 
discs. Each system evolves for 20 millions of years and the 
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initial conditions for each one are chosen random taking into 
account the following conditions: 

• The time-scale for the depletion of the gas, has a uni- 
form log distribution between 10 6 and 10 7 years. 

• The stellar mass has a uniform distribution in log scale 
in the range of 0.7 - 1.4M . 

• The distribution of metallicities of solar-like stars in 
the solar neighborhood follows a Gaussian distribution with 
fl = -0.02 and dispersion 0.22 (|Mordasini et alj|2009h . 

• The total mass of the disc is well approximated by a 
log-Gaussian distribution with mean —2.05 and dispersion 
0.85. We obtained this value by assuming a log-Gaussian dis- 
tribution and perfor med a non-linea r least squa re fit to the 
sampl e observed by lAndrews et al.l (|2009l ) and llsella et al.l 
(|2009h . 

• The characteristic radius, a c is also well approximated 
by a log-Gaussian distribution with jj, = 3.8 and a = 0.18. 
This distribution was obtained with the same procedure de- 
scribed in the previous item. 

We also note that the initial number of planets per disc 
depends on the disc size and mass as well as on the inner 
disc radius and the stellar mass, as is shown in Figure |3j 
where the initial number of embryos is plotted as a function 
of the disc mass for the three different values of 7 assumed 
in this work. Those discs characterized by a profile with 
7 = 1.5 are shown in figure 3(a) when 7 = 1 the results 



are shown in figure 3(b) and figure 3(c) represent the results 
when 7 = 0.5. 

We note in the figures, that the greater the mass of 
the disc, the lower the initial number of embryos. This is 
because the separation between the embryos is greater as 
the larger the initial mass of the embryo. On the other hand 
the initial mass of the embryo is larger when the density 
of solids is higher and this in turn, is higher as the greater 
the mass of the disc is. In summary, the greater the mass 
of the disc, the greater the separation between the embryos 
and consequently there will be a lower initial number of 
planetary cores. 

Figures [4] show the mass and semi major axis distribu- 
tion of all the planets that were formed in 1000 planetary 
systems, where the effect of planetary migration is not con- 
sidered and the initial density disc profile is 7 = 1.5 in 4(a) 



7 = 1 in 4(b) and 7 = 0.5 in the last figure. The first figure 



shows the results obtained with the profile more similar to 
the MMSN model, where the solids are available just beyond 
the snow line, fact that favours the formation of giant plan- 
ets at this location as we can see in the figure. This profile 
is also the steepest, making more feasible the formation of 
a single giant planet close to the snow line. In the second 
case, when 7=1 the profile is a bit soften, and as a conse- 
quence there are more solids available in the outer parts of 
the disc, fact that promotes the formation of giant planets 
further from the central star. Finally in the last case there 
are giant planets at around 15 au and even further. 

As a result it is noticed that in order to form giant plan- 
ets as the ones observed in our own Solar System in their 
current locations, it is more convenient to consider a proto- 
planetary nebula model with a softer profile. 

The next set of plots (figures [SJ) , are the results when 
we do consider planetary migration, but the type I planetary 
migration is delayed 100 times (c m i g i = 0.01). Figure 5(a 
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Figure 3. Initial number of embryos per disc as a function of the 
disc mass. The first figure (figure [3(c)[ l shows the initial number 
of embryos when the discs have a profile characterized by 7 = 1.5, 
when 7 = 1 the results are plotted in figure [3~(b)| and when 7 = 0.5 
the initial number of embryos per disc are shown in figure |3(c)| 



shows the distribution when the protoplanetary nebula was 
modeled with a density profile with 7 = 1.5, figure 5(b) 
represents all the planets formed when 7 = 1 and when 



7 = 0.5 we found the distribution observed in 5(c) 



Here we still note that the giant planet formation is 
favoured when a softer profile is considered. One conse- 
quence of the planetary migration is that a larger population 
of giant planets is observed. This is due to the fact that when 
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Figure 4. Mass and semi-major axis of all the planets formed 
in 1000 planetary systems, where the planetary migration was 
not considered. Figure |4(a)| shows the results when we consider 
7 = 1.5, figure [4(b) | shows the mass and location of the planets 
when 7 = 1 and finally the results obtained with 7 = 0.5 are 
shown in figure 14(c) I 



Figure 5. The figures show the distribution of mass and semi- 
major axis of all the planets formed when the planetary migration 
is considered but the type I is delayed 100 times. The first figure 
shows the planets formed when we consider 7 = 1.5, figure [5(b) | 
shows the results obtained with gamma = 1 and figurc [5(c)1 shows 
the distribution obtained with 7 = 0.5. 



a planet is migrating slowly it finds more solids available to 
accrete in the new feeding zone, and it grows faster. We also 
note that the more massive planets are found when 7 = 1.5, 
this is because this profile overestimates the mass of solids 
just beyond the snow line and the giant planets formed there 
have more solids available than in the other cases and as a 
result became bigger. Another point is that a population of 
giant planets closer to the star begins to be observed, fact 
that was not seen when the migration was not considered. 



This population of Hot-Jupiters will be called hereafter pop- 
ulation I. 

If the migration is delayed only 10 times (cmig = 0.1) we 



found the distribution shown in figure [U] Figure 6(a) shows 



the mass and semi major axis when 7 = 1.5, in figure 6(b) 
the profile is bit softer (7 = 1), and the results with 7 = 0.5 



are shown in figure 6(c 



Although migration is very fast, still can be seen that 
the smaller 7 is, the greater the chance of forming giant 
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region reached in solids grow faster and become giants in a 
time scale shorter than the type I migration time scale and 
starts to migrate with the slower type II migration regime, 
as a result, they remain close to the zone where they were 
originally formed. On the other hand we have the population 
of planets which where initially located inside the ice line, 
closer to the star, in this region the solids available are not 
enough to form giant planets faster than the type I time- 
scale, so they were moved rapidly inwards. On their path 
towards the star they found other embryos to collide with 
and new material to grow, so they became Hot-Jupiters. We 
note that when 7 = 0.5 the population is not well marked as 
in the other cases, this is because with this density profile, 
the giant planets are formed in the outer parts of the disc 
and there are no material available in the inner parts, so 
these discs do not favored the formation of hot Jupiters. 

Finally if we assume that there are no effects which act 
stopping or slowing down the migration rate (cmigi = 1), we 



observe the distribution shown in figures 7(a) (7 = 1.5), 7(b) 
(7 = 1) and 7(c) (7 = 0.5). In this cases the migration rate 



is faster than in the previous cases as a consequence some 
planets migrate really fast, faster than the depletion time 
scale, so the cores did not have enough time to grow and 
remain as small embryos. For this reason we note that the 
population of giant planets decreased in comparison with 
the cases of slower migration rate, and most of the planets 
are really close to the central star, located in the inner edge 
of the disc, specially when 7 = 1.5, because in this case the 
planets are formed closer to the star. 
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Figure 6. Mass and semi-major axis of all the planets formed in 
1000 planetary systems when the planetary migration was delayed 
10 times. Figure [6(a)] shows the results when we consider 7 = 1.5, 
figure |6(b)| shows the mass and location of the planets when 7 = 1 
and finally the results obtained with 7 = 0.5 are shown in figure 



planets farther away from the central star. In these distri- 
butions there are a more marked differentiation between two 
populations of giant planets that had begun to appear in fig- 
ures [S] the population I of giant planets close to the star, 
and another one of giant planets located beyond the ice line, 
farther of the host star. Since this differentiation is not ob- 
served in the distribution found without migration (figures 
|4}, this must be an effect caused by the orbital motion of 
the embryos. Those planets that were initially located in a 



3.1 Comparison with Exoplanets Observations 

In order to know which parameters must be considered to get 
a better representation of reality, in this section we compare 
the results obtained in our simulations with the mass and 
semi major axis distribution of observed exoplanets. 

Figure [5] shows the mass and semi major axis distribu- 
tion of observed exoplanets and the four giant planets in 
the Solar System, which are characterised for the biggest 
dots. We excluded of the sample those discovered planets 
that orbit binary and multiple stars, since the mechanisms 
of formation of these planets may not be explained with our 
model. 

We clearly notice in the Figure two populations: one is 
the population I and the another one are those exoplanets 
centered about ~ 3 au (population II, hereafter). 

Although in this work we are focused on the planetary 
migration as a possible explanation of these populations ob- 
served, it is not the only mechanism proposed when trying 
to explain how hot Jupiters came to be so close to their 
parent star: Kozai cycles and planet scattering should pro- 
vided us with an alternative explanation. These mechanisms 
excite the sky-projected obliquity of the planet and should 
provide us with a planet population on misaligned orbits 
with respect to their stars rotati on, that it is n o w beg inning 
to be detected. In the work of iTriaud et al.l (|2010h . they 
present observations of some exoplanets who present this 
misalignments. So this is saying that the population I could 
be explained by a combination of mechanisms. 



http: / /exoplanets. org 
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Figure 7. The figure shows the distribution of mass and semi- 
major axis found when we do not delayed planetary migration. 
When the density profiles corresponds to a value of 7 equal to 1.5 
we observe figurc [T(a)| figure [7(b) | corresponds to 7 = 1 and if we 
model the protoplanetary nebula with 7 = 0.5 we obtained the 
distribution observed in the last figure (Fig. |7(c)|. 



In order to get a better understanding of which are the 
parameters that gives us the best fit to the observations we 
superimposed the observations on our results found without 
considering the migration ^ and when it is delayed, con- 
sidering c = 0.01 (HO]), c = 0.1 ClTJ and c = 1 fl2]). for the 
three values of 7 analyzed in this work. 

In Figure [§] we show the overlap of the observations 
(exoplanets are grey dots and black big dots are Jupiter, 
Saturn, Uranus and Neptune) with the simulation results 



If 

w 




Semi major axis [AU] 

Figure 8. In the Figure the current mass and semi major axis dis- 
tribution of observed exoplanets is shown, where the giant planets 
in the Solar System are included(4 big dots) and the exoplanets 
found orbiting double or multiple stars were excluded of the sam- 
ple. 



when the effect of planetary migration is not considered. 
The case with 7 = 1.5 is shown in |9(a)| where the simula- 



1.5 is shown in 9(a) 
tion results fit pretty well the population II of exoplanets 
but it can not reproduce the formation of the giant plan- 
ets in our Solar System. When 7 = 1 there are more solids 
available to form giant planets farther from the star and 
as a result the population II of the observations can be re- 
produced and also the giant planets in the Solar System. 
Finally when 7 = 0.5 the simulations show that the giant 
planets are formed preferably about 10 au and although this 
allows the formation of the giant planets in the Solar Sys- 
tem, does not reproduce well the observed distribution of 
extrasolar planets. As a conclusion, the val ue 7 = 1, which 
is the approximate medium value found by lAndrews et al.l 
(2009) on their observations, are the best value in order to 
reproduce the population II of observed exoplanets and the 
giant planets in the Solar System, but it fails reproducing 
the population I of exoplanets. 

In order to try to form the planets belonging to pop- 
ulation I, we introduce the effects of orbital migration of 
planets. Figure \W\ shows the results for different values of 
7 when the migration is considered, but the type I migra- 
tion was slowed down 100 times (c m i g i = 0.01). Figure 10(a) 
shows the results with 7 = 1.5, when 7 = 1 we found figure 
and the third figure was obtained when considering 
0.5. We can notice that population I has began to ap- 



10(b) 



7 



pear, but this is stronger in the first (figure 10(a) I and in 



the second one (Figure 10(b) I 



In the next set of figures (111]) we show the overlap of the 
observations (plus the giant planets in the Solar System) and 
our simulation results when the migration is considered and 
the parameter for delaying type I migration is c m i g i =0.1. 
When 7 = 1.5 we observe figure [Tl(a)[ when 7 = 1 we ob- 
tained 11(b) and with 7 = 0.5 we found the overlap shown 
in figure 11(c) As was explained in the previous section, 
in these figures the same two populations noticed in the 
observations appeared, which are a consequence of plane- 
tary migration. Nevertheless these two populations are best 
differentiated in Figures [TT(a)1 and [TT(b)] which means for 
7 = 1.5 and 7 = 1. 
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Figure 9. In the figure we show the mass and semi major axis 
distribution obtained when planetary migration is not considered 
which are the small black dots, the distribution of observed exo- 
planets are the grey dots and the giant planets in the Solar System 
(the big black dots). Figure [9 (a) | shows the results of the simula- 
tions when 7 = 1.5, |9(b)"| when 7 = 1 and the results with 7)0.5 
are shown in|9(c)| 



Figure 10. Mass and semi major axis distribution of the planets 
found in the simulation when the migration is delayed 100 times 
are the small black dots, the current distribution of exoplants 
observed are the grey dots and the giant planets in the Solar 
System are the black big dots. Figure flO(a) | shows the results of 
the simulations when 7 = 1.5, 7 = 1 are shown in |10(b)] and the 
last figure shows the results with 7 = 0.5. 



Although these two figures show the results that are the 
best fit to the observations, the presence of a large popula- 
tion of very massive giant planets (with masses larger than 
5 M j U piter) is noticed here and it is not observed. If we cal- 
culate the percentage of observed exoplanets with masses 
greater than 5 Mj upite r and located at a distance < 0.4 au, 
we find that the 2.37% of the observed exoplanets have these 
characteristics, while the same percentage calculated for the 



results of our simulations shows that the number of super- 
hot Jupiters for the case where 7 = 1.5 and CMigi ~ 0.1 
is 30.86% and the percentage is 18.58% when 7=1 and 
Cmi g i — 0.1. We could think that these super hot- Jupiter 
population exists but it was not detected yet. But, if such 
planets exist, they should have been observed, since their ob- 
servation is favored by the observational techniques, so these 
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Figure 11. The figures show the mass and semi major axis of 
the giant planets in the Solar System, which are the biggest black 
dots, the distribution of observed extrasolar planets (grey dots), 
and the simulations results (the small black dots), when planetary 
migration is considered and delayed 10 times and 7 = 1.5 (11(a)) , 
7 = 1 (11(b)) and 7 = 0.5 (11(c)) . 



overpopulation must represent a limitation of our model. We 
will discuss this in section [3721 

If there is no factor delaying type I migration the over- 
lap of observations and simulation results are those shown 
in figure [121 The results obtained with 7 = 1.5 are shown 
in Figure 



in Figure 
7 = 0.5 



12(a) 



those obtained with 7 = 1 are represented 
12(b) and the last Figure show the results when 
In this case we note that the migration is too fast, 
and most of the planets reach the inner edge of the disc, so 



Figure 12. In the figures the mass and semi major axis distri- 
bution of our simulation results when the migration is considered 
and it was not delayed (c mig [ = 1) are shown with the small 
black dots, the distribution of observed exoplanets are overlapped 
in grey dots and also Jupiter, Saturn, Uranus and Neptune are 
shown with the biggest black dots. Figure [l 2 (a) | sh ows the simu- 
lation results obtained with 7 = 1.5, Figure 1 1 2 (b) | represents the 
results found with 7 = 1 and if 7 = 0.5 we obtained figure [T"2(c)| 



with this value of c m i g i we can reproduce the population I 
of exoplanets observed but the population II is absent. 

In summary, we find that when we do not consider plan- 
etary migration we can only reproduce population II and 
when migration is considered but it is not delayed, then 
only population I is reproduced. As a conclusion, planetary 
migration is a critical factor to explain the observed distri- 
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bution and we need to get a better understanding of the 
factors that act slowing down and stopping planetary mi- 
gration in order to get a better explanation of the planetary 
formation process. 

On the other hand one of the main objectives of our 
work is to show the importance of considering different pro- 
toplanetary nebula profiles on giant planet formation. If we 
analyze which is the disc model that better explains the ob- 
servations, we found that both 7 = 1.5 and 7=1 represent 
them quite well, but we choose 7 = 1 as the value that allows 
us to have an initial protoplanetary disc consistent with the 
observations and capable of reproduce the exoplanets obser- 
vational sample and also the four giant planets of our Solar 
System. 

3.2 About the overpopulation of massive giant 
planets close to the star 

As discussed in the previous section, although we obtained 



a good fit to the observations (Figure 11(b) I, we found 



an overabundance of planets with masses greater than 
5 Mj U piter and located very close to the central star, which 
are not detected observationally. In order to try to determine 
the cause of this overabundance of giant planets in the inner 
edge of the disc, we explored different alternatives. On the 
one hand we analyzed the possibility that such planets may 
have lost much of its atmosphere by erosion due to stellar 
winds and as a result most of it s primordial gaseous enve- 
lope has escaped. According to lLammer et alJ l|2003h . the 
mass loss of hot-Jupiters due to this effect is ~ 10 12 5 s" 1 
or ~ 1.5 x 10 -2 Mj up Gy~ , which is an important effect for 
objects with masses less or equal than Jupiter, but is irrel- 
evant for super hot-Jupiters and therefore does not explain 
the observed mass excess. 

On the other hand we also analyzed the possibility 
that the dissipation of gas in the nebula has a significant 
effect on these super hot-Jupiters. The gaseous disc dissi- 
pates between 1 and 10 million years l|Haisch et all l200ll ; 
iHillenbrandl 120051 ) due to different effects, one of them is 
the v iscous draining on to the central star (|Hartmann et all 
1 19981). although the photoevaporation is a lso very important 
|Hollenbach et allll99l ; IClarke et alj|200ll ). In this paper we 
considered a simple model of gas dissipation in the disc, 
where the gas dissipates exponentially but independently of 
the distance to the central star, while i n reality, the gaseous 
dissipation is much more complicated. IClarke et all (|200ll l 
have shown that a low level of photoevaporative mass loss 
from the disc in the region between 5 and 10 au can pro- 
mote the eventual rapid switch off of the disc inward of this 
radius, this coincides with the epoch at which the accre- 
tion rate through the disc falls as a consequence of viscous 
draining, to a level comparable with the photoevaporative 
mass-loss rate. As a result the inner disc empties faster than 
the rest of the disc. This phenomena could stop the type II 
migration of a giant planet already formed. So this could be 
one of the possible reasons of why we do not observe very 
massive giant planets close to that star. In order to explore 
this effect we need a more rigorous model for the disc evo- 
lution, which will be incorporated in our following works. 



Finally, when comparing figure 11(b) with other plots 
obtained with different values to delay the migration, there is 



no overpopulation when migration is slower (Figure 10(b) I, 



so this effect could also be an effect of the c m i g i value con- 
sidered. 

When considering simultaneously the standard core ac- 
cretion scenario for giant planet formation and gaseous mi- 
gration for the orbital evolution, we face a serious problem. 
For standard disc models, type I migration are primarily 
inward, and their timescale is at least 1 to 2 orders of mag- 
nitude shorter than the disc lifetime (10 6 — 10 7 yr), which 
means the cores would be accreted by the central star before 
they could build up any substantial gaseous envelope. There- 
fore, in order to avoid being swallowed by their central stars, 
the planet's migration must be significantly slowed down or 
stopped somehow, to this end we used the factor Cmigi, as 
was explained on section [2~3l But the question is, how much 
should we slow down planetary migration? 

A number of scenarios have been proposed to solve the 
non stop migration problem. iMenou fc Goodman! l|2004l ) ex- 
plore how the uncertainties in the structure of protoplan- 
etary discs may affect type I migration rates and found 
that it can be sig nificantly s l owed down at opacity tran- 
sitions in the disc, iLvra et all (|2009h show that as the sur- 
face density and temperature fall the planet orbital migra- 
tion and di sc depletion timescales eventuall y become com- 
parable, in IPaardekooper fc Mellemal (|2006T ) it was shown 
through three-dimensional, radiation-hydrodynamic simula- 
tions that planets could suddenly move outward as well as 
inward , depending on the local opacity, iMatsumura et all 
l|2007l ) proposed that the low-viscosity regions in protostel- 
lar discs (dead zones) ca n significantl y slow down planet 
orbita l motion, iTerqueml (|2003T ) and iNelson fc Papaloizoul 
l|2004h shown that including m agnetic fields may sl ow down 
or even stop type I migration, iMasset etT 

3 (|2006h showed 

that surface density jumps in the disc can trap low-mass 
protoplanets reducing the type I migration rate to the disc's 
accretion rate. We also add that irradiati on should lead to a 
hole formed by photo e vaporating wind l|Clarke et al.ll200ll : 
I Alexander et al.ll2006aH bh . which quickly depletes the disc, 
which possibly bring migration to an even earlier halt. 

All these works show that type I migration rate strongly 
depends on the exact disc conditions in the vicinity of the 
embedded embryo. Therefore we must be cautious about 
what is the c m igi value assumed. The factor that slows down 
migration should be a parameter that depends on the disc 
conditions in the proximity of the embryo and as such it is 
expected to be not a unique value for all discs but a distri- 
bution of possible values depending on the initial disc char- 
acteristics. We also propose that this could be a possible so- 
lution to the problem of the overabundance observed, since 
the number of massive giant planets observed in population 
I, is strongly dependent on the value of c m i g i considered, in 
other words, on the rate of type I migration assumed. So in 
order to have a more appropriate model, we should consider 
a distribution of values Cmigi, taking into account the fac- 
tors that inhibit the migration depends on the disc where 
the planets are embedded. 



4 SUMMARY AND CONCLUSIONS 

In the present work we have developed a semi-analytical 
model for computing planetary systems formation which is 
based on the core instability model for the gas accretion of 
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the embryos and the oligarchic growth regime for the accre- 
tion of the solid cores. The gas accretion model considered 
is a new approximation that we obtained introducing an 
analytic approxima tion to the numerical results found by 
l|Fortier et alj|2009h . 

The disc surface density profiles plays a fundamental 
role in the formation of giant planets and therefore, in the 
formation of the Planetary System, so in order to get a clear 
understanding of which are the initial conditions that al- 
lowed the formation of the observed exoplanets and the gi- 
ant planets in the Solar System we explore different models 
for the initial protoplanetary nebula structure. We assume 
that the gas and solid surface density are characterised by a 
power-law in the inner part of the disc and an exponential 
decay in the outer parts, which has great advantages, one of 
which is that we are not forced to arbitrarily cut the disc in 
order to obtained a finite disc mass. Both discs change with 
time: the solid disc is locally depleted due to the accretion 
of the embryos and the gas disc is globally depleted with a 
time-scale between 10 6 and 10 7 years. 

In our model we assume that the embryos have an or- 
bital motion due to their interaction with the host star and 
consider two regimes of planetary migration: type II for the 
larger embryos and type I for the smaller ones, where the 
rate is delayed a factor Cmi g i in order to represent the factors 
that act slowing down or even stopping planetary migration. 

With this model we consider different initial conditions 
to generate a variety of planetary systems and analyse the 
giant planets statistically. We assumed different initial discs 
profiles, a large range of discs masses and sizes according 
to the last results of observations in protoplanetary discs 
ijAndrews et alJ [20091 : llsella et alJl2009h . different stars in a 
range between 0.7 and 1.4Mq, and also different values of 
the parameter which acts delaying type I migration. We ex- 
plore the effects in the formation of assuming different discs 
and found that the observed population of exoplanets, in- 
cluding Jupiter, Saturn, Uranus and Neptune, is well repre- 
sented when considering planetary migration and a surface 
density with a power-law in the inner part characterised by 
an exponent of 7 = 1, which represents a softer profile when 
compared with the case most similar to the MMSN model 
case. Nevertheless we would like to note that despite having 
chosen 7 = 1 as the best value, in reality we should have 
considered a distribution of values of 7 in order to take into 
account the different discs where the planetary formation 
can occur. But, since the media value of the observations is 
7 ~ 1 (Andrews et al., 2009), the fact of not considering a 
distribution does not represent a big problem in the model 
considered. 

On the other hand if we explore the effects of consider- 
ing different values for delaying type I migration, we found 
that the best match to the observations is obtained when 
we delayed it some value between 10 and 100 times, but as 
noticed in section 13.21 this parameter depends strongly on 
the disc characteristics and therefore should be taken from a 
distribution of values and not be a unique value as has been 
considered so far in the synthetic population models. So we 
should have a better understanding of the factors that acts 
slowing down and stopping planetary migration in order to 
get a better explanation of the planetary formation process. 
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